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1. Introduction 

Although there are many good reasons to consider General Relativity (GR) as the 
best theory for the gravitational interaction, in the last few decades the advent 
of precision cosmology tests appears more and more to suggest that this theory 
may be incomplete. In fact, besides the well known problems of GR in explaining 
the astrophysical phenomenology (i.e. the galactic rotation curves and small scale 
structure formation), cosmological data indicates an underlying cosmic acceleration 
of the Universe which cannot be recast in the framework of GR without resorting 
to a additional exotic matter components. Several models have been proposed pQ 
in order to address this problem and currently the one which best fits all available 
observations (Supernovae la [2] , Cosmic Microwave Background anisotropies [3] , Large 
Scale Structure formation [A], baryon oscillations [5], weak lensing [6]), turns out to 
be the Concordance Model in which a tiny cosmological constant is present [7] and 
ordinary matter is dominated by a Cold Dark component. However, given that the 
A - CDM model is affected by significant fine-tuning problems related to the vacuum 
energy scale, it seems desirable to investigate other viable theoretical schemes. 

It is for these reasons that in recent years many attempts have been made to 
generalize standard Einstein gravity. Among these models the so-called Extended 
Theory of Gravitation (ETG) and, in particular, non-linear gravity theories or higher- 
order theories of gravity have provided quite interesting results on both cosmological 
[HI [HI El [H] and astrophysical [TUHH] scales. These models are based on gravitational 
actions which are non- linear in the Ricci curvature R and/or contain terms involving 
combinations of derivatives of R [13l HU [15] . The peculiarity of these models is related 
to the fact that field equations can be recast in such a way that the higher order 
corrections provide an energy - momentum tensor of geometrical origin describing an 
"effective" source term on the right hand side of the standard Einstein field equations 
[HI [9] . In this Curvature Quintessence scenario, the cosmic acceleration can be shown to 
result from such a new geometrical contribution to the cosmic energy density budget, 
due to higher order corrections to the HE Lagrangian. 

Among higher order theories of gravity, fourth order gravity f(R) models and 
in particular f(R) = R n [22] and f(R) = + [HI [HI [21] have recently gained 
particular attention since they seems to be able to provide an interesting alternative 
description of the cosmos [16]. Furthermore, these models can be related to other 
cosmologically viable models once the background dynamics has been introduced into 
the field equations [17j , providing a possible theoretical explanation to some of them. 

Because the field equations resulting from HTG are extremely complicated, 
the theory of dynamical systems provides a powerful scheme for investigating the 
physical behaviour of such theories (see for example J22[ [23]). In fact, studying 
cosmologies using the dynamical systems approach has the advantage of providing 
a relatively simple method for obtaining exact solutions (even if these only represent 
the asymptotic behavior) and obtain a (qualitative) description of the global dynamics 
of these models. Consequently, such an analysis allows for an efficient preliminary 
investigation of these theories, suggesting what kind of models deserve further 
investigation. 

In this paper, using the Dynamical Systems Approach (DSA) approach suggested 
by Collins and then by Ellis and Wainwright (see [24] for a wide class of cosmological 
models in the GR context), we develop a completely general scheme, which in principle 
allows one to analyze every fourth order gravity Lagrangian. Our study generalizes 
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[2"2"] . which considered a generic power law function of the Ricci scalar f(R) = R n and 
extends the general approach given in a recent paper [35] . Here a general analysis was 
obtained using a one -parameter description of any f(R) model, which unfortunately 
turns out to be somewhat misleading. 

The aim of this paper is to illustrate the general procedure for obtaining a phase 
space analysis for any analytical f(R) Lagrangian, which is regular enough to be well 
defined up to the third derivative in R. After a short preliminary discussion about 
fourth order gravity, we will discuss this general procedure, giving particular attention 
to clarifying the differences between our approach and the one worked out in [25j . In 
order to illustrate these differences and the problems that exist in |25j , we apply our 
method to two different families of Lagrangian R p exp qR and R + aR n . The last 
part of the paper is devoted to discussion and conclusions. 



2. Fourth Order Gravity Models 

From a purely theoretical point of view there are no prescriptions which prevent one 
to describe the gravitational interaction using a Lagrangian that is non-linear in the 
Ricci scalar and/or contains combinations of the Ricci and Riemann tensor. The 
main argument that led us to choose what we call the Hilbert - Einstein Lagrangian is 
that only in this case does one obtain second order field equations in the metric and 
the Newtonian Poisson equation in the low energy limit. If we relax the assumption 
of linearity in the gravitational action, the general coordinate invariance allows, in 
principle, infinitely many additive terms to the HE action: 

A G = J dPxy/^g [A + c R + Cl R 2 + c 2 R^R^ + c 3 R^ a pR^ aP + ....] , (1) 

and the general equations turn out to be particularly difficult to solve. However, if we 
limit ourselves to the fourth order, several simplifications are possible. First of all the 
Gauss-Bonnet theorem allows one to eliminate the R^v poR^ P<J terms. Moreover, in 
the case of homogeneous and isotropic spacetimes, the terms coming from the variation 
of the R I1U R' 11 ' invariant coincides with the one coming from the variation of the R 2 
term. Finally, one can define a suitable parametrization which makes it possible to 
recast the higher order field equations as a system of second order differential equations 
together with the constraint given by the definition of the new variables [35] . 

Consequently, in cosmology, the "effective" fourth order Lagrangian can be 
considered a generic analytic function of the Ricci scalar /(-R)ll: 

L = V~9 [f(R) + £m] • (2) 
By varying equation @, we obtain the fourth order field equations 

1 

T 



f(R)R^ - -f(R)g^ = f(Ry a0 {g^gpu ~ 9^9^) + f$ , (3) 



where f% = g_<5(^L M ) ^ ^ 

prime denotes the derivative with respect to R. 

^9 og^ 



Standard Einstein equations are immediately recovered if f(R) = R. When f'(R) ^ 
the equation ([3]) can be recast in the form 

1 

L fiu ^ flV ^ ^ flV ' 



GflU — RflV n9flvR — T^fiv ~ ^-flV "t" TflV ' (4) 



§ Unless otherwise specified, we will use natural units (h = c = fcs = 8nG = 1) and the (+. — , — , — ) 
signature. 
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where 

T *v = [l^u lf(R) - Rf'(R)} + f(Ry at3 bcvQfiu - > ( 5 ) 

represent the stress energy tensor of an effective fluid sometimes referred to as the 
"curvature fluid" and 

T M = 1 fM (6) 
fw f'(R) ' 

represents an effective stress-energy tensor associated with standard matter. 

The conservation properties of these effective fluids are given by the Bianchi 
identities Tj£. When applied to the total stress energy tensor, these identities reveal 
that if standard matter is conserved the total fluid is also conserved even though the 
curvature fluid may in general possess off-diagonal terms [221 123 12E] • In other words, 
no matter how complicated the effective stress energy tensor T^ OT is, it will always be 
divergence free if TM/ V = 0. When applied on the single effective tensors, the Bianchi 
identities read 

T$" f"{R) -j 

f'(R) f'(R) 2 



rpM;u _ _ J [_) T M f?iV (7\ 



f"(R) 

1 fiv — fUJAZ /»" 



the last expression being a consequence of total energy-momentum conservation. It 
follows that the individual effective fluids are not conserved but exchange energy and 
momentum. It is worth noting that even if the effective tensor associated with the 
matter is not conserved, standard matter still follows the usual conservation equations 

Let us now consider the Friedmann-Lemaitre-Robertson- Walker (FLRW) metric: 

dr 2 



ds 2 = dt 2 - a 2 (t) 



r 2 {d9 2 + sin 2 9d(j) 2 ) 



1 — kr 2 

For this metric the action the field equations |5]) reduce to 



(9) 



and 



* i n -i (10) 

2H + H 2 + = -- | - [f'R - f(R)] + f>- 3Hf> + p m j , 



R = -6 ( 2H 2 + H + A ) , (ID 



where H = a/a, f = and the "dot" is the derivative with respect to t. The 

system (fpQ|) is closed by the only non trivial Bianchi identity for T™ : 

(i m + 3if(/i m +p m ) = , (12) 

which corresponds to the energy conservation equation for standard matter. 
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3. The dynamical system approach in fourth order gravity theories 

Following early attempts (see for example [29]), the first extensive analysis of 
cosmologies based on fourth order gravity theory using the DSA as defined in [24] 
was given in [22] . Here the phase space of the power law model f(R) — xR 71 was 
investigated in great detail, exact solutions were found and their stability determined. 
Following this, several authors have applied a similar approach to other types of 
Lagrangians [30] , and very recently this scheme was generalized in [25j . 

In this paper we give a self consistent general technique that allows us to perform 
a dynamical system analysis of any analytic fourth order theory of gravity in the case 
of the FLRW spacetime. 

The first step in the implementation of the Dynamical System Approach is the 
definition of the variables. Following [22], we introduce the general dimensionless 
variables : 

/' R f 



" f'H' V 6H 2 ' ~ 6/'i? 2 ' ^ 

n = -£=L_ K = k 



3f>H 2 ' " a 2 H 2 ' 

where u m represents the energy density of a perfect fluid that might be present in the 
model [jj. 

The cosmological equations ([TO]) are equivalent to the autonomous system : 
dx 

e (2K + 2z - x 2 + {K + y + l)x) + fie (-3w - 1) + 2, 



dN 
dy 
dN 

dz 

dN 

dfi 

dN 
dK 

dN 



ye(2y + 2K + xq + A), 

ze (2K - x + 2y + A) + xeyq, (14) 
fle(2K - x + 2y-3w + 1), 
Ke(2K + 2y + 2), 



where N = \ lna\ is the logarithmic time and e = \H\/H. In addition, we have the 
constraint equation 

\ = -K -x-y + z + Sl, (15) 

which can be used to reduce the dimension of the system. If one chooses to eliminate 
K, the variable associated with the spatial curvature, we obtain 

= e (4z - 2x 2 + (z- 2)x - 2y) + fie (x - 3w + 1), 
^L= y£ [2n + 2(z + l)+x(q-2)}, (16) 

| In what follows we will consider only models containing a single fluid with a generic barotropic 
index. This might be problematic in treating the dust case because the condition w = might lead to 
additional fixed points. This issue has been checked in our calculations and no change in the number 
of fixed points has been found. In addition, the generalization to a multi-fluid case is trivial: one has 
just to add a new variable Q for each new type of fluid. This has the consequence of increasing the 
number of dynamical equations and therefore, the dimension of the phase space. However, since this 
generalization does not really add anything to the conceptual problem (at least in terms of a local 
analysis) . 
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(17) 



dz 

— = ze (2z + 2fl - 3x + 2)z + xe yq, 

^— = Q e (2tt - 3x + 2z - 3w - 1), 
K = z -\- VL — x — y — 1. 
The quantity q is defined, in analogy with [25], as 

a = JL 
Rf" ' 

The expression of q in terms of the dynamical variables is the key to closing the system 
([29]) and allows one to perform the analysis of the phase space. The crucial aspect to 
note here is that q is a function of R only so the problem of obtaining q = q(x, y, z, f2) 
is reduced to the problem of writing R = R(x, y, z, f2). This can be achieved by noting 
that the quantity 

r - -*f ( 18 , 
is a function of R only and can be written as 

r = - . (19) 

z 

Solving the above equation for R allows one to write R in terms of y and z and close 
the system (fTOj) . 

In this way, once a Lagrangian has been chosen, we can in principle write the 
dynamical system associated with it using (|16[) . substituting into it the appropriate 
form of q = q{y, z). This procedure does however require particular attention. For 
example, there are forms of the function / for which the inversion of (|19p is a highly 
non trivial (e.g. f(R) = cosh(i?)). In addition, the function q could have a non-trivial 
domain, admit divergences or may not be in the class C , which makes the analysis 
of the phase space a very delicate problem. Finally, the number m of equations of 
(|16p is always m > 3 and this implies that fourth order gravity models can admit 
chaotic behaviour. While this is not surprising, it makes the deduction of the non- 
local properties of the phase space a very difficult task. 

The solutions associated with the fixed points can be found by substituting the 
coordinates of the fixed points into the system 

H = aH 2 , a = -1-fli + Xi- Zi, (20) 
3(1 + w) 

Mm = T A*m , (21) 

a t 

where the subscript "z" stands for the value of a generic quantity in a fixed point. 
This means that for a^O the general solutions can be written as 

a =a (t-t ) 1/Q , (22) 

fi m = a (t - t ) ( « ' . (23) 

The expression above gives the solution for the scale factor and the evolution of the 
energy density for every fixed point in which a^O. When a = the (|2H)l reduces to 
H = which correspond to either a static or a de Sitter solution. 

The solutions obtained in this way have to be considered particular solutions of 
the cosmological equations in the same way that solutions can sometimes be found by 
using an ansatz for the form of the solution. For this reason it is important to stress 
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that only direct substitution of the results derived from this approach can ensure that 
the solution is physical (i.e. it satisfies the cosmological equations (fTU)) ). This check 
is also useful for understanding the nature of the solutions themselves e.g. to calculate 
the value of the integration constant (s). 

Also, the fact that different fixed points correspond to the same solutions is due 
to the fact that at the fixed points the different terms in the equation combine in such 
a way to obtain the same evolution of the scale factor. This means that although two 
solutions are the same, the physical mechanism that realizes them can be different 

In the following we will present a number of examples of f{R) theories that can 
be analyzed with this method and we compare the results obtained with those given 
in [25J [jj. An analysis of the approach presented in [25] and the differences with our 
method are given in the Appendix. 



4. Examples of f(R) - Lagrangians 

In this section we will show, with the help of some examples, how the DSA developed 
above can be applied. In particular we will consider the cases f(R) — R p exp(qR) 
and f(R) = R + aR n . Since the aim of the paper is to provide only the general 
setting with which to develop the dynamical system approach in the framework of 
fourth order gravity, we will not give a detailed analysis of these models. A series of 
future papers will be dedicated to this task. In what follows, we will limit ourselves 
to the finite fixed points, their stability and the solutions associated with them. A 
comparison with the results of [25] will also be presented. 



4.1. The f{R) = R p exp(qR) case 

Let us consider the Lagrangian f(R) — i? p exp(qi?). As explained in the previous 
section, the dynamical system equations for this Lagrangian can be obtained by 
calculating the form of the parameter q. We have 

q(y,z) = 2 V * . (24) 
y l — p z z 

% One difference betwen our approach ad the one in 1251 is that we consider a non zero spatial 
curvature k. This choice has been made with the aim of obtaining a completely general analysis of 
a fourth order cosmology from the dynamical systems point of view. In addition, since most of the 
observational values for the cosmological parameters are heavily model dependent, we chose to limit 
as much as possible the introduction of priors in the analysis. Anyway, the limit of fiat spacelike 
sections (K — > 0) can be obtained in a straightforward way for our examples. In fact, each fixed 
point is associated with a specific value of the variable K (i.e. a value for k) and the stability of these 
points is independent from the value of K. This means that in order to consider the limit K — » one 
has just to exclude the fixed points associated to K 7^ 0. Also, looking at the dynamical equations 
one realizes that K = is an invariant submanifold i.e. an orbit with initial condition K = will 
not escape the subspace K = and orbits with initial condition K 7^ can approach the hyperplane 
K = only asymptotically. As a consequence, one does not need to have any other information on 
the rest of the phase space to characterize the evolution of the orbits in the submanifold K = 0. 
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Substituting this function into (|2T)|) we obtain 
dx 



2x 



— = e [Az - 2x z + (z- 2)x - 2y] + ne{x-3w + 1), 

dy „ x z 

-j^=ye 2fl + 2z + 2 + - - 2 

dJSl y* — p z A 

dz 

dN y 2 — p 

dn 

dN 
K 



x y 



(25) 



Q e (2fl - 3x + 2z - 3w - 1) 
z + ft — x — y — 1. 



The most striking feature of this system is the fact that two of the equations have a 
singularity in the hypersurface y 2 = p z 1 . This, together with the existence of the 
invariant submanifolds y = and z — heavily constrains the dynamics of the system. 
In particular, it implies that no global attractor is present, so no general conclusion 
can be made on the behavior of the orbits without first providing information about 
the initial conditions. The finite fixed points can be obtained by setting the LHS of 
(|25[) to zero and solving for (x, y, z, f2), the results are shown in TableQ] 



Table 1. Fixed points of R p exp(<ji?).The superscript "*" represents a point 
corresponding to a double solution. 



Point Coordinates (x, y, z, ft) K 

~A (0,0,0,0) ~ 

B (-1,0,0,0) 

C (-l-3w,0,0,-l-3w) -1 

V (1 - 3w, 0, 0, 2 - 3w) 

£ (2,0,2,0) -1 

T* (1,-2,0,0) 

g (0,-2,-1,0) o 

H (4,0,5,0) 

J (2-2p,2p(l-p),2-2p,0) 2p(p-l)-l 

C* (-3(1 + w), -2,0, -4 -3w) 



u ( 4-2p (5-4p)p 5-4p n \ 

J 1 \ l-2p' 2p 2 -3p+l' (p-l)(2p-l) ' U J 

.f f -3(l+«;)(p-l) 3(l+tu)-4p -4p+3m+3 p(9w~2p(3w+4) + 13)-3(m + l) 



The solutions corresponding to these fixed points can be obtained by substituting 
the coordinates into the system ([20|) and are shown in Table [2][+|- 

The stability of the finite fixed points can be found using the Hartman-Grobman 
theorem [31] . The results are shown in Table [3] Note that some of the eigenvalues 
diverge for p = 0, 1. This happens because in the operations involved in the derivation 
of the stability terms p— 1 and/or p appear in the denominators. However this is not 
a real pathology of the method but rather a consequence of the fact that for these two 
values of the parameter the cosmological equations assume a special form. In fact, as 

+ Note that even if the parameter q is not present in the dynamical equations it appears in 
the solutions because we have calculated the integration constants via direct substitution in the 
cosmological equations. 
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Table 2. Solutions associated with the fixed points of W exp(gij). The solutions 
are physical only in the intervals of p mentioned in the last column. 



Point Scale Factor Energy Density Physical 



A 


a{t) = 


(* - to) 







V > 1 


B 


a(t) =a (t- t ) 1/2 







P > 2 


C 


a(t) = 


(t - t ) 







p > 1 


V 


a(t) =a (t- t ) 1/2 







P > 2 


e 


a(t) = 


(t - t Q ) 







p > 1 


IT* j 
J~ < 


a(t) = ao, 
a(t) = ao exp 


±^(t-to)_ 


n 
u 


P < 


p > o 

§,g > Vp > \,q < 


g 


a(t) = a , 
a(t) = ao exp 


±^#(*-*o) 





P < 


p > o 

|,g > Vp > |,g < 


n 


a(t) =a (t- to) 1/2 







p > 2 

i< p < i + yi 


i 


a(t) = (t- t ] 


^l-2p(p-l) 







c* 


a(t) = a , 
a(t) = ao exp 


_±^#(*"*o)_ 




) 


P < 


p > 

|,g > Vp > |,g < 


M 


a(t) = ao (t 


2p 2 -3p+l 


3(2p 2 -3p + l)(iu + l) 




n= I l s 

P 2 ' ' 4 


- to) 2 - p 


l-l>m f^mot p ~ 2 




N 


a{t)=a (t-io) 3(ro+1) 


Urn = f-moit — to) 2p 




- 4 (MmO 0) 



we will see in the next section, if one starts the calculations using these critical values 
of p one ends up with eigenvalues that present no divergence. 

Let us now compare our results with the ones in [25] (see the Appendix for details 
of this last method). The number of fixed points obtained for this Lagrangian, when 
K = 0, matches the ones obtained in This result can be explained by the fact 
that the solutions of the constraint equation for m (|A.1[) coincide with the ones coming 
from the correct constraint equation (|A.2[) (the matching between the two systems can 
be obtained setting w = in Table [T]). However, when one calculates the stability 
of these points our results are in strikingly different to those presented in [35] . For 
example, in our general formalism it turns out that the fixed point M is a saddle for 
any value of the parameter p and, as consequence, it can represent only a transient 
phase in the evolution of this class of models. Instead, in [25] the authors find that 
this point is a stable spiral and argue that this fact prevents the existence of cosmic 
histories in which a decelerated expansion is followed by an accelerated one. From 
this they also conclude that an entire subclass of these models (to = m(p) > 0) can 
be ruled out. Our results show clearly that this is not the case. Another example 
relates to the point Af. In [25] the authors find that this point can be stable when 
m — * 0, but from Tabled] one finds that this point can only be a saddle. As explained 
in the Appendix, the reason behind these differences is the fact that the method used 
in [25] leads to incorrect results when, like in this case, there is no unambiguous way 
of determining the parameter r = —y/z from the coordinates of the fixed points. 
Consequently the conclusions in [25] relating to the properties of these points are 
incorrect and have no physical meaning. 
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Table 3. The eigenvalues associated with the fixed points in R p exp(gR). The 
eigenvalues of the fixed point jV are displayed on three lines because of their 
mathematical complexity. 



Point Eigenvalues 



A 
B 
C 
V 
£ 
T 

Q 

n 

i 
c 

M 
N 



2, —2, 2, —1 — 3w] 

5,2,4,2-3w;] 

3(1 + -2, 2,l + 3«;] 

3(1 + 2,4,-2 + 3™] 
2,-2,2,-2-3(1 + 10 
4,-2,0,-4- Zw] 

•2,- 



1 / 68-25p 

2 V p-4 '' 



-5,4,2,-3(1 + ™)] 



68-25p 
p-4 ' ' 



-3(te + l) 



-2, -2 + p- v/3p(3p-4), -2 + p + ^3p(3p-4), 2p - 3(1 + to) 
-4,-2,0,4+3™] 



I + =±r, t^L, -2 + ^ + ^r, -4 



-1 - 



-1 



l-2p 1 p- 
p-3(i^+l) 
P 



l-2p 



. 6,2 ,1,3 

p-1 ' l-3p+p 2 ■ 
p-3{m + l) 
P 

3[(2p-l)tt)-l] _ j_ / -81(l+tu)+4p 2 (8+3^) 2 +3p(l+w)(139+87uQ-4p 2 (152+3^(55+18M>)) 

4p 4y (p-i)p 2 ' 

3[(2p-l)m-l] 1 / -81(l+m)+4p 2 (8+3tti) 2 +3p(l+u;)(139+87«;)-4p 2 (152+3m(55+18M))) 
4p + 4 V (p-l)p 2 _ 



& The f(R) = exp(qR) case 

Let us now consider the simple case of a pure exponential Lagrangian. Since this case 
has been extensively analyzed in a different paper [32j . we sketch here only the main 
results which are interesting for our discussion, referring the reader to [32] for further 
details. The function q is 

q(2/^) = £ , (26) 
V 

and dynamical system equations read : 
dx 

— = £ [4z - 2y - x(2 + 2x - z - fi)] + fie (1 - 3io) , 
J| = 2ye (2 + 2z + 2fi - z) , 

^ = 2zs(l + n + z-x), ( 2? ) 

— = Qe (2n - 1 - 3w + 2z - 3ar) , 

K =z + Q — x — y — 1. 

The coordinates of the fixed points, their eigenvalues and corresponding solutions are 
summarized in Table [5] . 

The cosmology of the Lagrangian f(R) = exp(qR) shows two interesting de Sitter 
phases: the point T> which is unstable and non hyperbolic point C that can behave as 
an attractor (see [32]). In addition it is possible to prove that there is a set of non zero 
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Table 4. The stability of the eigenvalues associated with the fixed points in the 
model R p cxp(qR). With the index + we have indicated the attractive nature of 
the spiral points. 

Point Stability 



A 

B 

C 

V 

£ 
T 

Q 

H 



C 
M 



saddle 

{repellor 
saddle 
saddle 

{repellor 
saddle 
saddle 
saddle 

attractor 
spiral + 
saddle 
saddle 

attractor < to< 1U i - ^ < p < V | < p < i + 
spiral + 
saddle 
non hyperbolic 

attractor < w< 1U p < ^(1 ,. ... . ., 

saddle otherwise 
saddle 



< w < 2/3 
otherwise 

2/3 < w < 1 
otherwise 



0<w<lU2<p<|| 
0<w<lU§<p<4 
otherwise 



V5 



0<w<lU0<p<§ 
otherwise 



V3)vi(l + V3) <p<2 



measure of initial conditions for which orbits connect these two points. In other words, 
such a Lagrangian can provide a natural framework both for inflation and the recent 
cosmic acceleration phenomenon. Nevertheless, it seems to lack an almost Friedmann 
phase which is required for structure formation. 



Table 5. Coordinates of the fixed points, the eigenvalues, and solutions for 
f{R) = exp(qij). The superscript "*" represents indicates a double point 



Point 


Coordinates [x, y, z, Cl) 


Eigenvalues 


Solution 




A 


[0,0,0,0] 


[-3ro-l,-2,2,2] 


a = a a (t — 


to) 


B 


[-1,0,0,0] 


[2-3^,2,4,4] 


a = a (t — 


tof* 


C* 


[1,-2,0,0] 


[-2,-4,-3™ -4,0] 


a = a e c ^~ 


-to) 


V 


[0,-2,-1,0] 


[-^±3,^,-2,-3-3^] 


a = a e c (*~ 


"to) 


£* 


[1 - 3w,0,0,2- 3w] 


[310-2,2,4,4] 


a = a a (t — 


tot 


T* 


[-3io-3,-2,0,-3«;-4] 


[3^ + 4,-2,-4,0] 


a = a Q e c (*~ 


-to) 


Q 


[-3w- 1,0,0, -3w- 1] 


[3w + l,-2,2,2] 


a = a a (t — 


to) 
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Table 6. The stability of the eigenvalues associated with the fixed points in the 
model exp(qR). See |32| for the stability of the non hyperbolic fixed points. 



Point Stability 



A 


saddle 


B 




f repellor < w < 2/3 


\ 


I saddle otherwise 


C 


non hyperbolic 


V 


saddle 


£ 




r repellor 2/3 < w < 1 




\ 


saddle otherwise 


T 


non hyperbolic 


Q 


saddle 



If we now compare our results with ones described in [25j . there are some clear 
differences. First, the number of fixed points turns out to be different. In fact, in our 
case there is no fixed point corresponding to the point P4 of |25j . and we obtain a 
new point T that does not appear in [25] . This follows directly from the pathological 
behaviour of equation (|A. 1|) . In fact, in the case of this Lagrangian, this expression 
and, in particular, the relation m(r) = — r — 1 has no solutions. Therefore, even in 
principle, there is no way to apply the method of [25] to this class of theories. Even if 
one refers to the correct equation (IA.2j) . the only possibility of having m(r) constant is 
to set x = 0, but this condition is not fulfilled by most of the fixed points. It is useful 
to see how these problems are related to the choice of taking m to be a parameter: 
if one substitutes the expression for m in terms of the dynamical variables into the 
initial dynamical system one obtains a set fixed points which correspond to the ones 
in Table El 

Finally, differences arise also in the stability analysis. The two points £ and T are 
non - hyperbolic and therefore require special treatment, without which it is impossible 
to draw general conclusions. Such treatment is given in detail in [32j . 

4.3. The case f(R) = R + aR n 

Let us discuss now the case of a Lagrangian corresponding to a power law correction 
of the Hilbert - Einstein gravity Lagrangian f(R) — R + aR n . In this case, the 
characteristic function q(r) reads: 
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and substituting this relation into the system of equations (I29[) one obtains 
= -2x 2 + {z- 2)x ~2y + Az + Q(x - 3w + 1), 

* =«(2*+2n-3* + 2)+ e -^, < 29 > 
dl\ n{z — y) 



fle(2n-3x + 2z-'3w- 1) 



dn 

dN 

K =z + ft — x — y— 1. 

As in the case of f(R) = R p exp(qR), the system is divergent on a hypersurface (this 
time y = z) but it admits only one invariant submanifold, namely y = and z = 0. 
This, again, implies that no global attractor is present and no general conclusion can 
be made on the behavior of the orbits without giving information about the initial 
conditions. The finite fixed points, their eigenvalues, their stability and the solutions 
corresponding to them are summarized in Tables [71 [SI [SI an d ITU1 



Table 7. Coordinate of the finite fixed points for R + aR n gravity. 



Point Coordinates (x, y, z, ft) K 

~A (0,0,0,0) ~ 

B (-1,0,0,0) 

C (-l-3w,0,0,-l-3w) -1 

V (1 - 3w, 0, 0, 2 - 3iu) 

£ (0,-2,-1,0) 

T (2,0,2,0) -1 

Q (4,0,5,0) 

H (2(l-n),2n(7i-l),2(l-n),0) 2n(n-l)-l 

T / 2(n-2) (5-4w)n 5-4n rA n 

X V 2n-l ' 2n 2 -3ri+l' 2rp -3n+l ' U J U 

r / 3(n-l)(tu+l) -4n+3m+3 -4ra+3m+3 -2(3m+4)n 2 +(9 W + 13)n-3(^+l) \ n 

*- I n ' 2n ' 2m 2 ' 2n^ / U 



As before our results are different from those given in 25J. First of all, out- 
set of fixed points do not coincide with the ones presented in [25 . In particular, in 
our analysis there is no fixed point corresponding to P§ a . Again, the reason for this 
difference is to be found in the constraint equation (|A.1[) , which in this case gives the 
incorrect set of solutions and therefore affects the set of fixed points. In fact, if one 
substitutes the expression for m(r) of [25] in terms of the coordinates in equations 
(34)-(39), it is easy to verify that two of these equations diverge at this point. 

The differences between the results in our approach and the one presented in |25j 
are even more evident when the stability analysis is considered. For example, the 
point £ , corresponding to Pi, is always a saddle, except into the region 32/25 < n < 2 
when it is attractive. This behavior is not obtained in [25] for which this point is 
stable only for —2 < n < 0. Also, points Q (P4) and T> (P3), which in our approach 
are always saddles (at least in the dust and in the radiation case), are always repellers 
in [55] . Finally, the stability of X corresponding to P§ appears to be different from the 
one presented in [25] . 
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Table 8. The eigenvalues associated with the fixed points in R + aR n 



Point Eigenvalues 



2,2,2,-3t0-l] 
5,4,2,2-3u;] 

2,2,3i0 + l,3(i0 + l) 
4,2,3(«; + l),3u) - 2] 

n 3 _ 1 / 25ra-32 
Z ' 2 2Y n >' 

-2,-2,2,-3( W +l)] 
-5,4,2,-3(1 + 10)] 

2(n 



25n-32 



-3(lB + l) 



1 ,n 



v /3n(3n-4), n - 2 + ^/3n(3n - 4), 2n - 3(w + 1) 



3[(2n-l)m-l] 



-81(l+tu)+4w 3 (8+3«;) 2 +3n(l+tu)(139+87M))-4n 2 (152+3«;(55+18u))) 
16(n-l)n 2 : 
3[(2n-l)w-l] , / -81(l+-m)+4n 3 (8+3u;) 2 +3K(l+tu)(139+87M))-4ri. 2 (152+3m(55+18m)) 
4n + V 16(n-l)n 2 ' 



3i0 



(5n+3(ri-2)tu-6) 



1 

n-1 i 



-l+2n ' 



i (3i0 



(5n+3(»-2)tu-6) 



3w 



5. Some remarks on the phase space of i?"-gravity 

In the previous sections we have analyzed two cases of fourth order gravity 
Lagrangians, and the relative subcases, to illustrate how a general dynamical system 
approach can be formulated for these theories. Furthermore, we have discussed the 
differences between our approach and the one presented in [25] . In this section we 
compare the results of these methods when they are applied to f{R) = xR n ■ The 
phase space of this class of theories has been investigated in detail in [22]. In the 
following we will show that only our method gives results that are in agreement with 
[22]. 

The crucial feature of i? ra -gravity in terms of the general method discussed above 
is that the characteristic functions r and q(r) are always constant. In particular, we 
have r = — n q(r) = n — 1. From the definition ([T^|) it is then clear that the variables 
z and y are not independent, i.e., the phase space i? ra -gravity is contained in the 
subspace y = nz of the general phase space described by (|29p . This can be easily seen 
if one substitutes y = nz into (|29p . Then the equations for y and z turn out to be 
exactly the same and (PI(?1) reduces to : 

= £ [-2x 2 + fix + zx - 2x - 2y + 4z + £1(1 - 3io)] , 

dy 
dN 
dQ 
dN 

with the constraint 



// ■ 2VL + \ ~2jx + 2(z + 1) 

Qe [2f2 - 3a; + 2z - 3w - 1] , 



(30) 



1 + x - y- z + K = , 



(31) 
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Table 9. The stability of the eigenvalues associated with the fixed points in the 
model R + aR n . The quantities A; related to the fixed point £, represent some 
non fractional numerical values (Ai fa 1.220, A\ fa 1.224, A 3 fa 1.470). 

Point Stability 



A 
B 
C 
V 



T 
Q 

n 

i 



c 



saddle 
repellor 
saddle 
saddle 
J repellor 
I saddle 
attractor 
spiral + 
saddle 
saddle 
saddle 
spiral + 
saddle 



< w < 2/3 
otherwise 

2/3 < w < 1 
otherwise 

i<™< 2 

0<n<§ 
otherwise 



< n < 1 
otherwise 




repellor 
saddle 

repeller 
saddle 



saddle, 
A± < n < A 2 U A 3 < n < 
otherwise 
1 < n < 



attractor 
repeller 
saddle 
repellor 
saddle 



n < 



1 < n < 



5 

4' 

otherwise, 

73) U n > 2, 

5 

4 ' 

otherwise, 
otherwise 



1(1 

1 < n < 



37), 



which is equivalent to the one given in [22] . Consequently the results obtained from 
the method presented above and [H] are identical Q 

The same cannot be said for the results of . In fact, although the set of 
fixed points are in agreement with the ones given in [22], the stability analysis is 
remarkably different. For example, the fixed point Q is claimed to become a stable 
spiral for n > 1, preventing the presence of orbits with transient almost-Fricdmann 
behavior. According to our results this is clearly not true since Q is always a saddle - 
focus or a saddle in such an interval of n. Furthermore, Q remains a saddle also 
for every n < 0.33 whereas in [25 is presented as a repeller for n — * — 1~. Other 
differences relate to the fixed point B, which is never stable in our approach, but is 
suggested to be attractive for 3/4 < n < 1 in [25] , 

* One can obtain a completely analogous result if the whole equations system (129 I t is considered 
without lowering the order of the equations. Of course one has to be careful in discarding the fixed 
points which do not fulfill the constraint y = nz. 
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Table 10. Solutions associated to the fixed points of R + aR n . The solutions 







JJIl_y SlCdl Ulliy 111 L11L- 111 LCI VeLlb <J1 fJ IIlCIlLlUIlCti 


111 LUC Ids Lr L-Ul U 11111. 




Point 




Scale Factor 


Energy Density 


Physical 


A 




a(t) = (t - t ) 





n > 1 


B 




a(t) = a (i - *o) 





n > 1 


C 




a(t) = (t - i ) 





n > 1 


V 




a(t) = a (t- t ) 1/2 





n > 1 




f o(t) = 


- ao, 




n > 


£* 


I = 


a exp [±2V3a 7 (2 - 3n)T(t - t )] , 





ra< §,a > V 






T 2(l-n) 




n > §, a < 


T 




a(t) = (t - to) 





n > 1 


Q 




a(t) = a (t - t ) 1/2 





n > 1 


H 




a(t) = y/l - 2n{n - 1) (t - t ) 





l<n>i + f 






2n 2 -3n+l 


3(2n 2 -3n+l)(TO+l) 




1* 




o(i) = fl (* - to) 2 ~" 


Mm = Mmot ™~ 2 


n = 5,Mm,o = 


C 




a(t) =a (t-t ) 3( " +1) 


Mm = Mm,o(t - to) 2P 


non physical 



Table 11. Coordinates of the fixed points for the model f(R) = X-R™. The 

superscript "*" represents a double solution. The point B is a double solution for 
n = 0,2. 

Point Coordinates (x,y, z,Q) K 

~A* (0,0,0,0) =T 

B (-1,0,0,0) 

C (-l-3ty, 0,0, -l-3u;) -1 

V (1 - 3to, 0, 0, 2 - 3w) 

£ (2(l-n),-2n(n-l),2(l-n),0) 2n(n - 1) - 1 

r I 3(w-l)(tu + l) -4n+3tu+3 -4rt+3m+3 -2(3m+4)n 2 +(9tti + 13)n-3(m+l) \ n 

I n ' 2n ' 2r^ ' 2r? ^ U 

2(n-2) (5-4»)n 5-4rt n \ n 

2n-l ' 2n 2 -3n+l> 2n 2 -3n+l > U J U 



Table 12. Coordinates of the correct fixed points for the model f(R) = xR n ■ 

Point Coordinates (x, y, z, Q) K 

~A (0,0,0,0) ~ 

B (-1,0,0,0) 

C (-l-3to, 0,0,-l-3u;) -1 

V (1 - 3to, 0, 0, 2 - 3w) 

£ (2(l-n),-2n(n-l),2(l-n),0) 2n(n - 1) - 1 

/ 3(»-l)(tu + l) -4n+3tu+3 -4rt+3m+3 -2(3m+4)n 2 +(9w + 13)n-3(u.+ l) \ „ 

I n ' 2n ' 27P ' §7? ^ U 

2(w-2j (5-4w)n 5-4n n \ n 

2n-l ' 2n 2 -3n+l' 2n 2 -3n+l ' U ,l U 
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Table 13. Stability of the fixed points for ii n -gravity with matter. We consider 
here only dust or radiation, see 22 for details and the case of stiff matter. The 
term "spiral+" has been used for pure attractive focus-nodes and the term "saddle- 
focus" for unstable focus-nodes. 



n<i(l- v / 3) i(l-V3)<n<0 < n < 1/2 1/2 < n < 1 



A saddle saddle saddle saddle 

B repulsive repulsive repulsive repulsive 

C saddle saddle saddle saddle 

V saddle saddle saddle saddle 

E saddle attractive spiral spiral 

T attractive saddle saddle attractive 



Kn< 5/4 5/4<n<4/3 4/3 < n < ±(1 + \/3) n>±(l + V3) 



A 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


repulsive 


repulsive 


repulsive 


C 


saddle 


saddle 


saddle 


saddle 


V 


saddle 


saddle 


saddle 


saddle 


H 


spiral 


spiral 


attractive 


saddle 


C 


repulsive 


saddle 


saddle 


attractive 



g n< 0.33 0.33<n<0.35 0.35 < n < 0.37 0.37 < n < 0.71 0.71 < n < 1 

w = saddle saddle-focus saddle-focus saddle-focus saddle 
w = 1/3 saddle saddle saddle- focus saddle- focus saddle- focus 

l<n< 1.220 1.220 < n < 1.223 1.223 < n < 1.224 1.224 < n < 1.28 

w = saddle-focus saddle-focus saddle-focus saddle-focus 

w = 1/3 saddle-focus saddle- focus saddle-focus saddle-focus 





1.28 < n < 1.32 


1.32 < n < 1.47 


1.47 <n< 1.50 


n > 1.50 


w = 


saddle-focus 


saddle 


saddle 


saddle 


w = 1/3 


saddle 


saddle 


saddle 


saddle 



6. Conclusions 



In this paper we have presented a general formalism that allows one to apply DSA to 
a generic fourth order Lagrangian. The crucial point of this method is to express the 
two characteristic functions [25 : 

q = w r = ~~ (32) 

in terms of the dynamical variables, which, in principle, allows one to obtain a closed 
autonomous system for any Lagrangian density f(R). 

The resulting general system admits many interesting features, but is very difficult 
to analyze without specifying the function q (i.e. the form of f(R)). Consequently, a 
"one parameter" approach can lead to a number of misleading results. 
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Even after substituting for q, the dynamical system analysis is still very delicate; 
in fact, q could be discontinuous, admit singularities or generate additional invariant 
submanifolds that influence deeply the stability of the fixed points as well as the global 
evolution of the orbits. 

After describing the method, we applied it to two classes of fourth order gravity 
models: R + aR n and R p exp(qi?), finding some very interesting preliminary results 
for the finite phase space. Both these models have fixed points with corresponding 
solutions that admit accelerated expansion and, consequently can be considered as 
possible candidates able to model either inflation or dark energy eras (or both). 
In addition, there are other fixed points which are linked to phases of decelerated 
expansion which can in principle allow for structure formation. These latter solutions 
are not physical for every value of their parameters, but this is not necessarily a 
problem. In fact, in order to obtain a Friedmann cosmology evolving towards a dark 
energy era, these points are required to be unstable i.e. cosmic histories coast past 
them for a period which depends on the initial conditions. This means that the 
general integral of the cosmological equations corresponding to such an orbit will only 
approximate the fixed point solution and this approximate behavior might still allow 
structures to form. 

It is also important to mention the fact that even if one has the desired fixed 
points and desired stability, this does not necessarily imply that there is an orbit 
connecting them. This is due to the presence of singular and invariant submanifolds 
that effectively divide the phase space into independent sectors. Of course one can 
implement further constraints on the parameters in order to have all the interesting 
points in a single connected sector, but this is still not sufficient to guarantee that an 
orbit would connect them. The situation is made worse by the fact that, since the 
phase space is of dimension higher than three, chaotic behavior can also occur. It is 
clear then, that any statement on the global behavior of the orbits is only reliable if 
an accurate numerical analysis is performed. However, these issues (and others) will 
be investigated in more detail in a series of forthcoming papers. 

A final comment is needed regarding the differences between our results and the 
ones given in [23]. Even if the introduction of q and r, was suggested for the first 
time in that paper, the results above (and in particular the existence of a viable 
matter era) are in disagreement with the ones given in that paper. The reason is that 
the authors of [25] used "a one parameter description" in order to deal with ([29]) in 
general. We were able to prove that, unfortunately, not only are the equations given 
in [23] incomplete, but also that the method also gives both incorrect and misleading 
conclusions. 

Appendix A. The approach of [25j 

The basic idea for closing the general system of autonomous equations for /(iZ)-gravity 
was suggested for the first time in [25] . In fact, if we define m(r) = q _1 the equations 
([29]) for w = and K = are equivalent to the ones given in this paper. The authors 
of [25] proposed that the function m could be used as a parameter associated with the 
choice of f(R), thus obtaining a "one parameter approach" to the dynamical systems 
analysis of f(R) gravity. Unfortunately their method has several problems that lead to 
incorrect results. These problems can be avoided only if one considers the framework 
presented above. 
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Let us look at this issue in more detail 0. In [25] the system equivalent to 
associated with the relation 

dr R 

1N = r(l + m(r) + r)—, (A.l) 

which is clearly a combination of the equations for z and y. In order to ensure that 
the variable r and consequently the parameter m is constant they require the RHS 
of the above equation to be zero. Their solution to this problem is the condition 
1 + m(r) + r = 0, which is an equation for r when the function m(r) has been 
substituted for and is also the bases of their method of analysis. 

The problem here is that this equation has not been fully expressed in terms of 
the dynamical system variables. In fact, one can rewrite (|A.1[) in the form: 

dr r(l + m(r) + r) 



dN m(r) 
at the condition 
r(l + m(r) + r) 



-x , (A.2) 



dv 

which means that the condition = in fact corresponds to 

dN 



x = (A.3) 
m (r) 

rather than 1 + m(r) + r = 0. Equation (|A.3[) has a solution if 
x = 0, 

r = °' (A.4) 
(1 + m(r) +r) _ 

m(r) 

and this leads to solutions for r which are in general different from the values of r 
obtained from 1 + m(r) + r = 0. This inconsistency has major consequences for the 
rest of the analysis in [25] , leading to changes in the number of fixed points as well as 
their stability (see the text above for details). 

In fact, a more careful analysis reveals that for some of the fixed points (e.g. 
Pi,.. .Pi) the values of r obtained from the relation r = —y/z either cannot be 
determined unambiguously or do not solve the condition 1 + m(r) + r = 0, which 
is claimed to come from (|A.1[> in [25] . 

This is a clear indication that the approach used in [25] is both incomplete and 
leads to wrong conclusions. It is also interesting to stress that if one substitutes the 
expression for m in terms of the dynamical system variables in (26-29) of [25 , the 
results match the one obtained in our formalism. This implies that the reason the 
method described in [25] fails has its roots in the attempt to describe the phase space 
of a whole class of fourth order theories of gravity with only one parameter. 

(j It is important to note that in 1251 the signature is not the same of the one used here (e.g -,+,+,+ 
instead of +,-,-,-) and the definition of the variables are slightly different. The transformation from 
one variable to another is as follows: 

x — > — asi, y — > —X'i, z — > X2, K — > 0, w — > 0. 

However, as expected, this does not affect our conclusions. 
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